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In evaluating neoclassical transport by radially-local simulations, the magnetic drift tangential to a flux surface 
is usually ignored in order to keep the phase-space volume conservation. In this paper, effect of the tangential 
magnetic drift on the local neoclassical transport are investigated. To retain the effect of the tangential 
magnetic drift in the local treatment of neoclassical transport, a new local formulation for the drift kinetic 
simulation is developed. The compressibility of the phase-space volume caused by the tangential magnetic 
drift is regarded as a source term for the drift kinetic equation, which is solved by using a two-weight Sf 
Monte Carlo method for non-Hamiltonian system [G. Hu and J. A. Krommes, Phys. Plasmas 1, 863 (1994)]. 
It is demonstrated that the effect of the drift is negligible for the neoclassical transport in tokamaks. In 
non-axisymmetric systems, however, the tangential magnetic drift substantially changes the dependence of 
the neoclassical transport on the radial electric field E r . The peaked behavior of the neoclassical radial fluxes 
around E r = 0 observed in conventional local neoclassical transport simulations is removed by taking the 
tangential magnetic drift into account. 


I. INTRODUCTION 

Neoclassical transport caused by Coulomb collisions in torus plasma is fundamental for a magnetically confined 
plasma since it determines an irreducible minimum for the plasma transport. It also plays a key role in determining the 
radial electric field through the ambipolar condition of the neoclassical particle flux when non-axisymmetric devices 
such as stellarators and heliotrons are considered. In addition, the neoclassical viscosity caused by non-uniform 
magnetic field influences plasma parallel flows. 

The neoclassical transport theory is based on the drift kinetic equation, in which the fast gyration of the plasma 
particle is removed. Many analytic and numerical evaluations have been done for axisymmetric tokamaks and non- 
axisymmetric devices. 1-8 For this purpose, additional assumptions are usually made in the drift kinetic equation. At 
first, the higher order radial drift is neglected. This enables ones to solve the “radially local” drift kinetic equation, 
leading to “local” neoclassical transport, where “local” means that the drift kinetic equation and the neoclassical 
transport is only determined by its radially local parameters. Second, the tangential component of the magnetic drift: 

v B = v B - (v B ■ Vs) e s (1) 

is omitted, where «b is the magnetic drift composed of the VS drift and the curvature drift, s is a label of magnetic 
flux surfaces, and e s is the covariant basis vector in s direction. Third, the mono-energetic particle assumption is of 
importance. This assumes that the particle velocity v , or kinetic energy mv 2 /2 is unchanged along the particle orbit. 
Finally, E x B drift is assumed to be incompressible to conserve the phase-space volume. With these assumptions, 
the evaluation of the neoclassical transport becomes much easier since the drift kinetic equation described in five¬ 
dimensional phase space is reduced to that in three-dimensional phase space. 

As mentioned above, the neoclassical transport simulations are based on many assumptions, which are interde¬ 
pendent. The main purpose of this paper is to reconsider the validity of the approximations, especially with respect 
to the tangential magnetic drift, v B . Depending on the approximations made in the drift kinetic equation, various 
kinds of the drift kinetic equation and the particle orbit appear in this paper. (A) The drift kinetic equation without 
all the assumptions described above. Since the equation includes the higher order radial drift term in this case, the 
neoclassical transport with the finite orbit width (FOW) effect can be evaluated. 9 The neoclassical transport is also 
called a global one due to the fact that it involves the radially global effect in it. It should be noticed that since v is 
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in proportion to the product of the radial drift and the radial electric field E r , we also call v the FOW effect in this 
paper. (B) The drift kinetic equation without the radial drift term. Since the radial drift term is neglected, the drift 
kinetic equation becomes local, and the local neoclassical transport is obtained. We refer to this particle orbit as zero 
orbit width (ZOW) orbit in order to distinguish it from other kinds of the local orbit. (C) The drift kinetic equation in 
ZOW limit with vb = 0. The particle orbit and the drift kinetic equation in this limit have many preferable features, 
as described later in Sec. II, some authors evaluate this type of the local neoclassical transport. We would like to call 
it the zero magnetic drift (ZMD) limit. (D) ZMD limit with the mono-energetic particles and incompressible E x B 
drift. In this limit, the particle orbit reduces to the same one that is adopted in a widely-used neoclassical transport 
code, DKES. 7,10 We call this particle orbit as the DKES-likc orbit. 

Conventionally, the neoclassical transport has been evaluated locally, and DKES-like orbit has been adopted in many 
codes. This is justified in a typical torus plasma if the radial drift term is negligibly small. Recently, however, several 
authors have pointed out that there are some cases where the approximations in the conventional local neoclassical 
transport models are violated. For example, the FOW effect becomes significant near the axis of a tokamak due to 
the potato orbit. 11 The electron FOW affects the neoclassical transport in a high electron temperature stellarator due 
to its complicated orbit and low collisionality. 12 Also, the mono-energy assumption may cause underestimation of the 
fraction of the helically-trapped particles in a quasi-symmetric stellarator when the radial electric field is finite. 13 

So far, efforts have been made to investigate the effects of the FOW and/or the mono-energetic particle, while the 
effect of the tangential magnetic field on the local neoclassical transport has not been considered . 12,14 " 15 Although 
conventional local neoclassical codes provide reliable results in many cases, there seems to be several situations when 
the magnetic drift needs to be included, e.g., a resonant behavior of the magnetic drift with E x B drift. The 
resonant behavior called the poloidal resonance between the tangential magnetic drift and E x B drift occurs in a 
non-axisymmetric magnetic field configuration. Dependence of the neoclassical transport on the radial electric field 
is qualitatively varied by the poloidal resonance. However, the influence of the tangential magnetic drift on the 
neoclassical transport is not fully clarified since there are no local neoclassical transport models which include the 
effect. This makes it difficult to compare neoclassical transport models with and without the tangential magnetic drift. 
For example, when comparing the global neoclassical transport to the local one in DKES-likc limit, the compressible 
E x B drift, finite v and Db in addition to the effect of the radial motion simultaneously affect the neoclassical 
transport of the global model. In other words, there exists a large gap between the global neoclassical transport 
model, in which the effect of the tangential magnetic drift is included, and conventional local models ignoring the 
effect. It is necessary to bridge the gap by developing a local neoclassical transport model based on the drift kinetic 
equation in the ZOW limit in order to explore the effect. 

In this paper, we present a new formulation of the local drift kinetic simulation in the ZOW limit, where the 
tangential magnetic drift term is retained while the radial drift term is ignored. Due to the tangential magnetic 
drift, the phase-space volume, and thus the particle number are not conserved. The resultant local drift kinetic 
equation becomes non-Hamiltonian, and the compressibility of the phase-space volume acts as a source term. Hu 
and Krommes prescribes the two-weight Sf method appropriate for such non-Hamiltonian system with an arbitrary 
source/sink term. 1 ' 1 Based on their work, we develop a numerical code for the local neoclassical transport with vb- 
The code developed here requires less computational cost than global ones with the FOW effect due to the local 
feature of the code. It provides a more accurate method to evaluate the neoclassical transport in plasmas where the 
tangential magnetic drift becomes significant. Another advantage of the code is that the particle orbit in the code 
can be switched among the ZOW, ZMD and DKES-like limit models. This enables us to investigate the effect of the 
particle orbit on the local neoclassical transport. We investigate the neoclassical transport in an axisymmetric and a 
non-axisymmetric plasmas using the code. It is found that the neoclassical transport in the ZOW limit is almost the 
same as that in DKES-like limit in an axisymmetric case. This suggests that vb is not significant in axisymmetric 
tokamak as expected. On the other hand, the local neoclassical transport in ZOW limit is demonstrated to show 
a radial electric field dependence, or the poloidal resonance at a finite E r , which has not been seen in the local 
neoclassical transport. When using the ZMD and/or DKES-like limit orbits, a large peak of the neoclassical radial 
flux is observed at E T = 0. Such a large neoclassical transport is removed in the ZOW limit due to the effects of the 
tangential magnetic drift. As a result, the neoclassical transport in the ZOW limit approaches to that in the global 
FOW model. 

The remaining part of this paper is organized as follows. The drift kinetic equations with FOW effect and in 
various local limits (ZOW, ZMD, and DKES-like limits) are described in Sec. II. The property of the phase-space 
conservation in each limit is also presented. The two-weight 5f Monte Carlo method for non-Hamiltonian system is 
described in Sec. III. Numerical results for axisymmetric case and non-axisymmetric case are presented in Sec. IV 
and V. A summary is given in Sec. VI. 
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TABLE I. Comparison of guiding-center orbits and conservation properties included in the global and local drift kinetic models. 
Model equations for models (A), (B), (C) and (D) are described in corresponding subsections of Sec. II. In the table, Comp. 
and Incomp denote compressible and incompressible, respectively. 



Global 

Local 

Model 

(A) FOW 

(B) ZOW 

(C) ZMD 

ZMD 

(D) DKES-like 


(Finite Orbit Width) 

(Zero Orbit Width) 

(Zero Magnetic Drift) 

+ mono-energy 


Particle orbit 

Full orbit 

finite Db 

D B = 0 

vb = 0, 
mono-energy 

v B = 0, 

mono-energy, 

Incomp. E x B 

A 

0 

finite 

0 

finite 

finite 

V 

oc —e&ip 

oc —e&ip 

oc —e&ifi 

0 

0 

vb 

Included 

Included 

None 

None 

None 

V ■ z 

0 

finite 

0 

finite 

0 

E x B drift 

Comp. 

Comp. 

Comp. 

Comp. 

Incomp. 

Dimensions 

5 

4 

4 

3 

3 


II. DRIFT KINETIC EQUATION FOR NEOCLASSICAL TRANSPORT 


We derive several kinds of the drift kinetic equation for the first order distribution function /i based on various 
models of the guiding-center particle orbit stepwisely in the following subsections. Our starting point is the drift 
kinetic equation without any assumptions, which is adopted in the the radially-global neoclassical transport with the 
FOW effect, such as GTC-NEO 17 and FORTEC-3D 9 codes. The radially local drift kinetic equation is obtained by 
omitting the higher order radial drift term, sdfi/ds, from the equation, leading to the ZOW orbit. Then, further 
simplifications are usually made to the local drift kinetic equation in conventional neoclassical transport simulations 
instead of solving the ZOW-limit equation directly. In ZMD limit, a term involving the tangential magnetic drift to a 
flux surface is approximated to be zero, that is, Vb • V/i = 0. In addition, mono energetic particle vd v f\ oc sd s f\ = 0, 
and incompressible Ex B drift are assumed in DKES-like orbit, where v is the particle velocity. Subsidiary changes are 
also introduced in some essential conservation properties of the drift kinetic equations along with these assumptions 
for the local neoclassical transport. The differences among these global (FOW) and local neoclassical transport models 
are summarized in Table I for the later convenience. 


A. Drift kinetic equation with finite orbit width (FOW) effect 


Staring equation is the drift kinetic equation for the guiding-center distribution function, f a = f a (R,v,f): 18 


dfa ; • d f a 

¥^SI- C(U 


( 2 ) 


where subscript a represents the particle species and the phase-space variables are denoted by 2 as = (R, y,£); R 
is the position vector, and £ = v\\/v is the pitch angle of the parallel velocity U|| = v ■ b with the unit vector parallel 
to the magnetic field, b = B/B-, C{f a ) is the linearized collision operator acting on f a . The subscript a is omitted for 
simplicity hereafter unless it is necessary. In the following, we use Boozer coordinates 19 to specify the position vector 
R as R = (if, 0,(), where is the toroidal magnetic flux, 9 and ( are the poloidal and toroidal angle variables, the 
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drift equations of motion are derived from the canonical Hamiltonian given by White : 20 
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(3b) 

(3c) 
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(3e) 


where m and e are the mass and electric charge of the species; t{ip) is the rotational transform; $ = $>{ip) is the 
electrostatic potential; prime denotes the derivative with respect to ip; 7 = G{ 1 + + /(« — yy^G") is used for 

simplicity with the poloidal and toroidal current fluxes, G{ip) and I {ip). 

By introducing a small parameter S ~ 0{p/L ), where p represents the Larmor radius and L denotes the typical 
scale length, the drift kinetic equation can be solved order by order. To this end two important orderings are assumed; 
one is the transport ordering of ~ 0{5 2 u>t), where to t ~ L/vth is the transit frequency, and the other is the drift 
ordering of fE/^th ~ 0(5), where ue and ve is the E x B drift and its magnitude, and Vth is the thermal speed of 
the particle. For the drift ordering, we use the fact that ve is tangential to a flux surface since $ is a flux function. 
By decomposing the distribution function as / = /o + / 1 , it can be readily shown that /o is Maxwellian, that is, 
/o = fu{ip,v) from the leading order drift kinetic equation. 

The drift kinetic equation for /1 becomes as follows: 


Eh 

Dt 


dfi j 5/i 
~dt + ^W 
j df 0 . df 0 


)dfi :3/i - d ^ 1 cdfi 

l ar + c ar + ' , 8^ + ? ar 


( 4 ) 


where /o = fo{ip, v) is used, and C(f 1 ) = C(fi, fo) + C{fo, / 1 ) is a linearized collision operator for Coulomb collisions. 
It should be noted that eq.(4) involves terms of different orders, ipd^fi and vd v f\ are of the order of O{6 2 ujtfo), 
where v oc 4 Yip. On the other hand, Odgfi, (d(fi and are composed of O{6u> t fo) term arising from the parallel 

motion and O{5 2 ojtfo) terms from the perpendicular drift. It should be noted that d t fi is regarded as the order 
of O{6 3 ojtfo ) hr the quasi-steady state according to the transport ordering. Solving eq. (4) directly with linearized 
collision operator for /1 leads to the neoclassical transport with the FOW effect which is represented by the radial 
drift term ip and velocity term v. 

The drift kinetic equation with the higher order radial drift, satisfies the conservative properties of phase-space 
volume and particle number. This is due to the fact that the guiding-center motion in the five-dimensional phase- 
space, z, is Hamiltonian. For i, Liouville’s theorem is satisfied: 


■z = ^J2{jzi) = 0 , 


J=i 


( 5 ) 


where z J = {ip, 9, £, v, £), and the Jacobian of the five-dimensional phase-space coordinates is given as 


2 nB*v 2 G + tI 
B B 2 ’ 


( 6 ) 


with £?j| = ^B + V x bj ■ b. Here, we assume that the Jacobian does not depend on time explicitly. Therefore, 

the original (global) drift kinetic equation, eq. (2) and (4), conserves the five-dimensional phase-space volume. 

The particle number in the phase-space is also conserved. This can be readily seen by rearranging the five¬ 
dimensional drift kinetic equation in the conservative form. The operator V is defined as 


V = 


dt 


1 _ E. 

Jdzt J dt 


- 


dzi 


( 7 ) 
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where Liouville’s theorem, eq. (5), and dj/dt = 0 are used to show the second equality. Using V and assuming 
collisionless limit of G(/) —► 0, the drift kinetic equations for / and f\ becomes 


Vf = 0 


Vfi 


} df Q . df 0 


Integrating the drift kinetic equations over the entire phase-space volume dz 
above, we have 


( 8 ) 

(9) 

Jdz 1 • • • dz 5 and using the definitions 


dN 
dt 
dN i 
dt 


= 0 


= 0 . 


( 10 ) 

( 11 ) 


In the equations above, the total particle numbers for / and f\ in the phase-space, are defined by N = f dzf and 
Ni = J dzfi , respectively. The particle number conservation also holds for collisional cases since a proper choice of a 
collision operator satisfies the conservation laws for the particle number, momentum and energy . 9 


B. Local drift kinetic equation in the Zero Orbit Width (ZOW) limit 


The local drift kinetic equation in ZOW limit is obtained by neglecting the higher order radial drift term, ifdfi/dif, 
in eq. (4). We have 


Dfi 

Dt 


dfi xdh l3/i <9/i -<9/i 

-W +e ~M +< lK +v ~^ +( ^( 

j df 0 . df 0 . 


( 12 ) 


It should be noticed that D/Dt represents the total derivative along the particle orbit in z^ = (0,£,v,£). not in 
z in eq. (4). What is important in this equation is that the radial variable if only appears in the right hand side 
of the equation as a source term, ipdfo/dif. This means that the dependence on if only enters through ifdfo/dif as 
a parameter when equilibrium distribution f Q is given. The radially local neoclassical transport at a surface can be 
evaluated independently by solving eq. (12) at the surface. It should be noted that, in ZOW limit, the higher order 
radial drift term, ifd s fi, is ignored while vd v f\ term still remains in order to make the comparison to the further 
reduced local neoclassical transport model (the ZMD limit) simpler. 

The particle orbit lies on a single flux surface during the time evolution. The drift equations of guiding-center 
motion for the local drift kinetic equation are obtained by omitting the effect of the radial drift in eqs. (3) as 


6 
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(13a) 

(13b) 

(13c) 

(13d) 


To derive eqs. (13), B*^ ~ B is used, and the coefficient 7 is approximated as G + tl. Our main purpose of this paper 
is to construct a proper numerical method to solve eq. (12) along the particle orbit given by eqs. (13). 

The consequence of the neglect of the radial drift in the local drift kinetic equation, eq. (12), is compressibility of 
the phase-space volume. The divergence of the phase-space flow becomes finite due to the presence of the magnetic 
drift in poloidal and toroidal directions, vb- In contrast to the case of eq. (4), the conservative form V in z^ does 
not agree to the total derivative along the particle orbit in z^ -coordinates, D/Dt ; V = D/Dt + V • z/^. In this 
phase-space, the Jacobian is written as follows: 
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^B^v 2 G + tl 
B B 2 


~ 27 TV 2 


G + tl 
B 2 ’ 
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(14) 
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where B*^ is again approximated by B and the (G + 1 1) /B 2 part is Jacobian of Boozer coordinates and it is the same 
as that of the five-dimensional phase-space. It should be noted that this approximation of B jj to B does not influence 
on the conservation property of the phase-space volume, although the Hamiltonian nature of the system is broken. 
Using eq. (13), the compressibility of the phase-space volume can be obtained as 


_ mv 2 (1 + e) f 3 dB f dB dB\ 

2eB(G + tI) \B dip V d( d9 ) 

( d 2 B d 2 B\\ 

+ \ G diPd9 dipdQ J J ' 


(15) 


The right hand side of eq. (15) exclusively arises from ub part of 9 and f. Since (V • i (4 ^) /i is of the order of 0(S 3 ), 
it is a higher order contribution to the local drift kinetic equation, (12). The local drift kinetic equation, (12), is 
rewritten as follows: 


Vh = (V • i (4) ) h + S 0 (16) 

where Sq = — 1S use( i to formally represent its behavior as a source term in the local drift kinetic equation. 

Integrating eq. (16) over the phase space gives rise to the non-vanishing contribution to dNi/dt: 

d -£=Jc iz< 4 >(v. *<*>)/,, (17) 

where contributions from other source term, So vanishes since they represent the radial velocity moment of Maxwellian 
distribution function, fo . The number of particle in the phase-space is not conserved in the local drift kinetic equation 
due to the compressibility when the finite tangential magnetic drift is considered. The same situation occurs for the 
conservation of the magnetic moment /q the magnetic drift contributions in 9 and C, again lead to the violation of 
p = 0 : 


M = 


mv±v± 
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mv 2 , ( ■ dB ■ dB \ 
2 B 2 \ 6 ~d9 + ^~dC ) 


B dip ^ 


(18) 


Although the phase-space volume, and thus the particle number Ni are not conserved when considering the ZOW- 
limit particle orbit, this does not cause any matter practically in evaluating steady-state neoclassical transport ob¬ 
servables such as the particle and energy fluxes in many cases by the Sf Monte Carlo method prescribed in Sec. III. 
In fact, the neoclassical transport observables presented in Secs. IV and V reach steady-state values in our particle 
simulations while Ni remains negligible compared to N. The ZOW model is inappropriate only when extremely large 
radial excursion of the guiding centers, such as the potato orbit near the axis, 11,21 mainly determines the neoclassical 
transport. In fact, all the local models of the neoclassical transport are insufficient in such cases, and the global, or 
the neoclassical transport with the FOW effect is essentially required. It should be also noted that, although V • 
acts as a source term in the local drift kinetic equation, the term is quite different from the source term of Landreman 
et al. 15 It is pointed out in the reference that, for the cases of the drift kinetic equation with the full and partial 
trajectories, surface averaged conservation laws of the particle number and energy result in a singular perturbation 
problem when E r approaches to 0, where the full trajectory corresponds to the ZMD orbit in this paper. Their source 
term is introduced to remove the singular perturbation. On the other hand, V • term is introduced here to solve 
the local drift kinetic equation with the non-Hamiltonian property by particle simulations. 

The radial locality and the requirement of the divergence-free phase-space flow do not hold simultaneously. It 
should be noticed that this cannot be avoided even if one chooses other variables for the velocity space. For example, 
when we chooses (u|| , q) as independent variables instead of (u,£) and keep p = 0, then the conservation of the total 
energy along the particle orbit does not hold. Nevertheless, as presented later in Sec. Ill, one can numerically solve 
the local drift kinetic equation by treating the compressibility as another source term in a system as well as So- 


C. Local drift kinetic equation in the Zero Magnetic Drift (ZMD) limit 

In some neoclassical transport simulations, additional assumption is made for the local drift kinetic equation in 
ZOW limit presented in the previous subsection; Db ■ V/i =0 is assumed to neglect the tangential magnetic drift. 
Among them are EUTERPE code 22 and the works of Landreman, 14,15 for example. Hence we refer to this type of 
the particle orbit as a ZMD orbit to distinguish it from the ZOW orbit and DKES-like orbit in this paper. 
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By omitting the magnetic drift from eq. (13), the equations of the guiding-center drift are given as 
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(19a) 

(19b) 

(19c) 

(19d) 


The local drift kinetic equation is formally the same as eq. (12). Solving eq. (12) along with this modified guiding- 
center orbit, eqs. (19), leads to the local neoclassical transport without Db- 

V • = 0 is again recovered due to the absence of the magnetic drift terms in 9 and £ and the presence of v term 

in this limit. (Neglecting v term again gives rise to the finite V • z^\ see eq. (21).) This leads to the conservation 
of the particle number N\ since conservative form V agrees with the total derivative along the orbit, D/Dt. The 
conservation of /j, is also followed from the absence of Db- It should be noted that, although the magnetic drift terms 
in 9 and £ are of the order of S 2 as well as the E x B drift, the latter is only taken into account in this limit. This 
causes a large peak of neoclassical radial fluxes around E r = 0. 


D. Mono-energetic assumption and incompressible E x B drift (DKES-like limit) 


In this subsection, the local drift kinetic equation in the ZMD limit is further reduced by assuming the so-called 
mono-energetic particles (v = 0), in which the kinetic energy does not experience any change. As shown later in this 
subsection, the assumption of the mono-energetic particle again violates the phase-space volume conservation. To 
recover the conservation property, the E x B drift is also assumed to be incompressible, resulting in the DKES-like 
limit particle orbit. It should be noted that we always use the DKES-like limit (v = 0 and incompressible E x B 
drift) when considering the mono-energetic particle assumption in this paper. 

When the mono-energetic assumption, vdfi/dv = 0, is made in addition to the assumption of Db ■ V/i = 0, the 
local drift kinetic equation reduces to three-dimensional problem, Since v term is higher order as described above, 
this assumption is consistent to neglecting the higher order radial drift. Under the mono-energetic assumption, the 
local drift kinetic equation (12) becomes 


d Jl + e d J .i + + 

dt do ^ d c 


= S 0 + C(h). 


( 20 ) 


In the equation, the particle velocity v (kinetic energy mv 2 / 2) only enters parametrically through dfo/dv in the 
right hand side. Also, the test-particle collision operator Ct(/i) included in the linearized collision operator C(/i) 
should be also modified under the mono-energetic assumption. The test-particle collision operator is reduced to the 
pitch-angle scattering operator (Lorentz operator). The particles do not experience the energy scattering. The four 
dimensional phase space z^ reduces to three-dimensional one, z^ = (9,(,0- v can be treated just as a parameter 
to solve the equation as well as the radial variable s. The mono-energetic guiding-center drift equations of motion is 
the same as eqs. (19) except for v = 0 in this case. 

The mono-energetic assumption again violates phase-space volume and particle number conservations. This arises 
due to the presence of the E x B drift in 9, £ and £. The divergence of z becomes 


v B 3 ) _ 1 3(1+e 2 ) d$ 

J 2 B 3 dil> 




( 21 ) 


where J is Jacobian given in eq. (14). E x B drift also results in fi oc <Af> /dip\ the magnetic moment is not conserved 
along the mono-energetic guiding-center orbit. 

E x B drift is often regarded as an incompressible drift in conventional local neoclassical transport codes, such as 
DKES. 7,10 The E x B drift, is approximated as follows: 
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or equivalently, <J>' in the guiding-center drift equations of motion is replaced by & B 2 / (B 2 ). With this replacement, 
the guiding-center drift equations of motion then becomes 
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where Db = 0 and v = 0 are also assumed. While the local drift kinetic equation along this guiding-center orbit still 
violates the phase-space volume conservation, it is satisfied if effect of <f>' is simultaneously removed from £■ Indeed, 
this is what DKES and many conventional neoclassical transport code assume in their approach; no magnetic drift, 
mono-energetic particle, incompressible E x B drift, and no <f>' effect on £. We call this particle orbit DKES-like 
orbit, for simplicity. 

Adopting these all assumptions, the drift equations of motion for DKES-like orbit are 


e =ohi{ G W)' il+ ' v " B ) 

i= chi{- I W) ¥+v,B ) 

■ w (l — £ 2 ) fdB 0B\ 
~ ~ 2 {G +tl) V5C +t 06 ) 


(24a) 

(24b) 

(24c) 


As a consequence, the phase-space volume conservation is again satisfied with V • = 0. Hereafter, in this paper, 

the guiding-center particle orbit with incompressible E x B drift represents those described by eq. (24), not eq. (23). 
On the other hand, however, the violation of /i = 0 is not recovered even if incompressible E x B drift is assumed; p 
is not conserved along the guiding-center trajectory due to the radial electric field. 

These additional assumptions are simultaneously adopted in many conventional neoclassical transport codes. This 
makes it difficult to properly compare the difference among various models of the local drift kinetic equation and/or the 
drift kinetic equation with FOW effect. In order to address the effect of each drift on the local neoclassical transport, 
it is necessary to construct a numerical method to solve the wide varieties of the local drift kinetic equations with 
several drifts included/neglected independently. 


III. 2-WEIGHT Sf MONTE CARLO METHOD FOR LOCAL NEOCLASSICAL TRANSPORT 


Two-weight Sf Monte Carlo method is widely used to solve the drift kinetic equation and its formulation for 
collisional transport with incompressible flow of V • z = 0 was given in detail by Brunner et al. 23 and Wang et al 24 
respectively. Since our interest is the local drift kinetic equation in ZOW limit, where V • zf'O ^ 0, the formulation 
needs to be modified to appropriately treat such case. 

Hu and Krommes pointed out in their work, 16 the two-weight Sf Monte Carlo method is applicable to a non- 
Hamiltonian system in which the compressibility of the phase-space volume is included as a source term in the weight 
evolutions. According to the work, we apply the method to the local drift kinetic equation in ZOW limit. Below, we 
briefly review the standard formulation of the two-weight Sf Monte Carlo method for Hamiltonian (incompressible 
flow) system. The formulation is given for five-dimensional phase-space coordinates for generality. Then, to discuss 
the ZOW limit in four-dimensional case, the effect of V ■ zf'O term is included as a source term. Cases of ZMD and 
DKES-like limit are then presented. 

In the two-weight method, two weights, w and p, are assigned to each simulation marker. Then, the discretized 
distribution function of simulation markers, F = F(z, w,p\ t), is introduced. It is noted that the distribution function 
F is defined not in an ordinary phase-space z, but in an extended phase-space ( z,w,p ). Using F, the distribution 
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functions, /o and /i, are evaluated weighted sum of F as follows: 



F = ^25(z - Zi)S(w - Wi)5{p - pi)J x (z) 

i 

(25a) 


g= f Fdwdp = ^ S(z — z/)J~ x (z ) 

i 

(25b) 


fo= f pFdwdp = y ^Pi6(z - z i )J~ 1 (z) 

i 

(25c) 


fi = f wFdwdp = ^2 WiS(z - z-,)J a (z), 

i 

(25d) 

where g = g{z) denotes the marker distribution function in the ordinary phase space z, 
marker indices. The expressions for the weights Wi and pi are obtained by integrating eqs 
and g: 

and subscript i denotes the 
. (25) for f 0 and fi using F 


fl(Zi) 

Wi = 

g{Zi) 

fo(Zi) 

Pi= , , • 
g{ z i) 

(26a) 

(26b) 


The total derivative along the particle orbit including the test-particle collision D^/Dt is defined by rewriting 
eq. (4) as 


D (c) fi 

Dt 


Dfi 

Dt 


Ctp(/i) 


Sq + Cfp(/m), 


(27) 


where the linearized collision operator C(/i) is decomposed into the test-particle part, Ctp(/i) and field-particle 
one, Cfp(/m)- Since the simulation markers are discretized, the total derivative along the particle D^/Dt should 
be replaced by what is appropriate for the discretized markers. In the two-weight Sf Monte Carlo method, this is 
enabled by approximating the test-particle collision in D^/Dt by Monte Carlo collision operator for the discretized 
markers; 25 ’ 26 D ( m )/Dt ~ D^/Dt. It is worth noting that this approximation of the collision operator is the origin 
of the weight spreading. 23 

The marker distribution function g(z) is conserved along D^/Dt: 


D (c) g 

Dt 


Vg-g(V-z)-C TP (g)=0, 


since V • z = 0. Thus, we obtain following equation from eq. (27), 


D (c) fi 

Dt 


D (c >< 
’ Dt 


£>(e) 


w 


£>(c) 


W 


Dt 


= g- 


Dt 


(28) 


(29) 


Using D(°)/Dt ~ D/Dt for the right hand side of the second equality, the time evolution of w along the marker 
orbit with the Monte Carlo collision is obtained as 


D( M h 

Dt 


= --p- f V'TTT + V-7T - Cfp) /m 


d 


/m V d' 1 !’ 9v 


(30) 


where eq. (27) is used for the right hand side. Similarly, time evolution of p can be read as 


= + (31) 

When we consider the local drift kinetic equation in ZOW limit in z^ phase space, the compressibility of the phase- 
space volume remains in eq. (28). Hereafter in this section, we use formally the same notations for the distribution 
functions /i, g, etc., although they are defined in z^ phase space not in z. According to this change, the total 
derivatives D^/Dt and /Dt also become those defined in z^. 
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With finite V • z^\ the marker distribution along D^ /Dt becomes 


D {c) g 

Dt 



(32) 


For eq. (12), a similar discussion as in eqs. (27) - (29) leads to the time evolution of w as follows: 

* = + - Cfp ) /u+ ” ( v ■ i<4> ) (33) 

To obtain the time evolution of p, it should be noticed that the total derivative D^/Dt ~ /Dt is described in 
2 ;( 4 ) = (0, y, £). For p we obtain 


. 1 £> (c) /m 

Pi = -n -r~ 

g Dt 

_ Pi . df M 

/m V dv 


+ p(V-z (4) ) 

+ p(v-i (4) ) . 


(34) 


Thus, only v derivative and compressibility appear in the right hand side. The solution of the local drift kinetic 
equation (12) is obtained by following the time evolution along the orbit defined by eqs. (13) with the Monte Carlo 
test-particle collision. 

The source term of the phase-space incompressibility is of the order of /i (V ■ i^ 4 )) ~ w (V • - O(S 2 uj t f 0 ). 

This means that such a non-conservative property introduced to the local drift kinetic equation induces higher order 
effect on the neoclassical transport as the E x B drift, mono-energetic particle assumption, etc. Thus, the use of the 
finite Db can be justified in solving the local drift kinetic equation, eq. (12), which is of the order of O(6ojtfo). The 
violation of the conservation of the phase-space volume occurs in the local neoclassical transport calculations, if we 
would like to treat the finite v b in the local drift kinetic equation (ZOW limit). 

To avoid the phase-space volume compressibility in the z^\ v-q = 0 must be assumed. The phase-space volume is 
conserved along the orbit, that is, V ■ z/ 4 ^ = 0. The second terms in the right hand side of eq. (33) and (34) reduce 
to zero in the ZMD limit. 

Finally, in DKES-like limit, same arguments are made in three-dimensional phase-space coordinates, z^ = ( 6 , £, £), 
due to the mono-energy assumption. Noting that V • z^ = 0 in this phase space, time evolutions of the weights are 
described as 


* =_ l)('^ +< ’lr CFp ) /u (35) 

Pi = 0, (36) 

where the time evolution of p can derived by using the fact that the total derivative D <c ' ) /Dt is described in z^. The 
second weight pi of each simulation marker is conserved during a simulation. The two-weight Sf method in DKES-like 
limit can be interpreted as the one-weight method due to the mono-energetic particle. 


IV. AXISYMMETRIC CASE 

As a code verification, several neoclassical transport values are compared to theoretical estimations for an axisym- 
metric configuration. To see the difference among the particle orbit discussed in Sec. II, we use three kinds of the 
particle orbit and collision; one is DKES-like orbit with the pitch angle scattering and without the field-particle colli¬ 
sion operator (denoted as DKES-like, PAS); the second one also has the DKES-like particle orbit with full test-particle 
collision operator and field-particle collision operator (DKES-like, FC); the particle orbit of the third one is ZOW 
orbit with full test- and field-particle collision operators (ZOW, FC). For the reference against the radially-global 
neoclassical code, numerical simulations are also performed by using the global code, FORTEC-3D (F3D). It should 
be noted that the same full collision operator including the field particle operator as that described above is used in 
the global code FORTEC-3D. 

We use an axisymmetric tokamak geometry with a circular cross-section, of which equilibrium is constructed by 
VMEC code 27 with parameters below; Rq and a are the magnetic axis and minor radius are given as Rq = 2.35 m and 
a = 0.47 m, respectively; the aspect ratio at the plasma edge e~ 4 = R 0 /a = 5.0; the magnitude of the magnetic field 
at the magnetic axis is B 0 = 1.91 T. The safety factor q = 0.854 + 2.184p 2 is used. Hereafter we use the normalized 
toroidal magnetic flux, sj^/ip a , as a flux-surface label, where is the toroidal magnetic flux at the plasma edge. 
The plasma density n and ion temperature T | are shown in Fig. 1. The normalized collisionality, +(, is shown in 
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FIG. 1. The radial profiles of and ion temperature (Ti) and the plasma density (n e ). 



v/v a 


FIG. 2. The radial profile of the banana-normalized collisionality i/£ and the safety factor q. 

s 


Fig. 2, where is the normalized collisionality defined v £ = qR/{y th,i£ 3/,2 Tii); Uth,i = \fiT\Jrn\ denotes the thermal 
velocity of the ion; t;; represents the ion-ion collision time 3 defined from the ion collision time of Braginskii 28 t\ as 
Ti = \[2t \i In the Fig. 2 (b), the safety factor q is also represented. The radial electric field is set to be constant during 
a simulation and is given as a numerical parameter according to the force balance relation: 3 


<V5,|| B) 


gGT • 

e 


(dtP\ 

1 1 |A NC - 

^ d In Ti 

d In rii 

05 

a- 
I ^ 

\dr) 

i (Pi 

2 dr 

dr 

Ti dr J 


(37) 


where the ion parallel flow Vf| = 0 is assumed to be zero, and the coefficient /3; NC is also given in eqs.(6.134) and 
(6.135) in the reference. 3 It should be noted that /3; NC depends on the collisionality; since the collisionality is artificially 
varied to see the collisionality dependence in numerical results presented below, the radial electric field is also varied 
depending on the collisionality. The E T determined as such enables us to obtain the steady-state neoclassical transport 
with less computational time. This does not affect numerical results as the neoclassical transport in an axisymmetric 
tokamak is independent of E T . 

The radial profile of the neoclassical ion thermal diffusivity yi is compared to theoretical estimates from Chang- 
Hinton formula 29 and the moment method of Hirshman-Sigmar 2,30 in Fig. 3. It should be noticed that the theoretical 
estimations for ion thermal diffusivity are obtained by imposing a certain limitation on the aspect ratio on the local 
drift kinetic equation along with the mono-energetic particles and zero tangential magnetic drift. It is demonstrated 
that the numerical results of ZOW and DKES-like orbit with full collision operator (FC) well reproduce the theoretical 
results over the wide region of the plasma, while the DKES-like orbit only with pitch angle scattering (without field- 
particle part) tends to underestimate the thermal conductivity. The result shows that the momentum-conservation of 
the collision operator more influences neoclassical transport simulations than the difference in the particle orbits. On 
the other hand, vb hardly influences on the neoclassical thermal transport. This is because trapped particles with 
Db in the axisymmetric tokamak just precess in the symmetry direction, and this causes no additional neoclassical 
transport. 

The results of the DKES-like orbit with both collision cases (FC and PAS) increase towards the magnetic axis of 
y/ip/ip a < 0.15 as Chang-Hinton and Hirshman-Sigmar theories predict. On the other hand, the result of the ZOW 
orbit show decreases towards the axis after showing an unphysical increase around there. This suggests that the 
incompressibility of the phase-space volume in the ZOW limit, which is caused by the radial deviation of the particle, 
becomes significant. In fact, Xi °f the FOW case (denoted as global in the figure) also shows a discrepancy from those 
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FIG. 3. The radial profile of the ion thermal diffusivity obtained by simulations with several kinds of the particle orbit and 
collision. DKES-like denotes the particle orbit in the DKES-like limit, and ZOW is that in the Zero Orbit Width limit, and 
Global represents the radially global orbit obtained by FORTEC-3D. PAS represents that it uses only the pitch angle collision 
and does not include field particle operator; FC represents that the simulation uses the full collision operator composed of the 
test particle collision operator with the pitch angle and energy scattering and the field particle operator. Theoretical estimations 
from Chang-Hinton (C-H) formula and the moment method of Hirshman-Sigmar (H-S) are also shown by dotted and chain 
lines, respectively. 



FIG. 4. The radial profile of /3 ; NC obtained by simulations with several kinds of the particle orbit. Abbreviation for the kinds 
of the orbit and collision is the same in Fig. 3. Theoretical estimation of Hirshman-Sigmar (H-S) formula is also shown by a 
chain line. 


of the DKES-like orbit cases and theoretical values near the axis, where Xi smoothly decreases towards zero. The 
decreasing tendency towards the axis can be attributed to the existence of the potato orbit which has the large radial 
deviation near the axis of a tokamak. 11,21 The width of the potato o rbit be comes ~ 0.13 for the parameters used here, 
showing that the FOW effect due to the orbit becomes significant i/i/’/V’a. < 0.13. In other words, the compressibility 
of the phase-space volume introduced by the radial drift results in the unphysical transport there. It should be noted 
that Xi °f the global case shows a decrease in the edge region since the particle loss at the last closed flux surface is 
only included in the radially global simulation. Also, both theoretical predictions underestimate towards the edge 
region due to the effect of the finite aspect ratio, which is only partly included in the theoretical calculation. 

The radial profile of /3; NC is shown in Fig. 4. In the figure, /3; NC , is compared to theoretical values of Hirshman- 
Sigmar. 2 It should be noticed that the momentum conservation does not hold for the result of the DKES-like PAS 
case due to the absence of the field particle operator is not included in the simulations. The cases with full collision 
operator with both local orbits (DKES-like and ZOW) again show a better agreement with theoretical values. 

Then, the collisionality dependence of x'i and /3; NC at ^/^/'0 a — 0.49 is investigated. For this purpose, the collision- 
ality in Fig. 2 is numerically magnified by 0.01, 0.5, 5 and 10. It should be noted that E r given for each collisionality 
case is varied due to the difference of /3; NC at the initial state as described before. The results of Xi and /Sf 10 are 
shown in Figs. 5 and 6, respectively. It is shown that the numerical results of the ZOW and DKES-like orbit cases 
with the full collision operator (FC) show better agreement with both theoretical values over the wide range of the 
collisionality. of the ZOW and DKES-like orbits reproduce Hirshman-Sigmar estimates, especially in the low 
collisionality regime of < 0.1. 
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v b 


FIG. 5. Collisionality dependence of \i at p ~ 0.49. The same abbreviation as in Fig. 3 is used for the kinds of the orbit and 
collision. Chang-Hinton (C-H) and Hirshman-Sigmar (H-S) estimations are represented by dotted and chain lines, respectively. 
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FIG. 6. Collisionality dependence of /3; NC at p ~ 0.49. The same abbreviation as in Fig. 3 is used for the kinds of the orbit and 
collision. A chain line represents the Hirshman-Sigmar (H-S) formula. 


V. NON-AXISYMMETRIC CASE 

To see the effect of the tangential magnetic field drift, the neoclassical transport in an non-axisymmetric magnetic 
field configuration is investigated. For this purpose, we take an LHD configuration as an example. Due to the 
asymmetry in the magnetic field, the intrinsic ambipolar condition is broken; the neoclassical transport depends on 
the radial electric field. The dependence of the neoclassical transport on E I shows a resonant peak at a finite E T 
when the tangential magnetic drift exists, while the peak appears at E T = 0 without the tangential drift. 12 However, 
this was demonstrated by comparing local codes (GSRAKE 8 ’ 31 and DCOM/NNW 32 ) in DKES limit and a global 
code (FORTEC-3D), in which the FOW effect and the tangential drift was both included. The role of the tangential 
magnetic drift in the local drift kinetic equation is numerically studied below by using the local code developed here 
with the ZOW, ZMD and DKES-like orbits. The validity of the numerical results are also checked by comparing the 
results to conventional local neoclassical codes and FORTEC-3D. 

The so-called inward-shifted magnetic field configuration of LHD is chosen as magnetic axis R ax = 3.6 m, and 
the magnetic field strength at the axis is B ax = 3.0 T. The ion temperature T) = 1.0 keV and the density n e = 
0.5 x 10 19 m -3 at the axis. The equilibrium magnetic field is again constructed by a widely-used equilibrium code for 
three-dimensional field, VMEC. The plasma collisionality and the rotational transform are shown in Fig. 7. 

The E t dependence of the neoclassical particle flux at ^/ip/ip a ~ 0.29, 0.49 and 0.74 are shown in Fig. 8 - 10. E r is 
given as a constant parameter for each simulation. The local orbit in ZOW (w/ i>b), ZMD (w/o ub) and DKES-like 
limits are used in evaluating the flux by the local code developed in this paper. The full collision operator including 
the field-particle one is used for all orbit types. The particle flux by DKES, GSRAKE and FORTEC-3D (denoted 
as Global, F3D) are also shown in the figures. The former two are the local codes, while the latter is the global 
code. DKES code used here includes the momentum correction. 33 GSRAKE solves the bounce-averaged drift kinetic 
equation with the local DKES-like orbit, and only the pitch angle scattering collision operator without the momentum 
correction is used. As mentioned above, FORTEC-3D code evaluates the neoclassical transport with the FOW effect 
and Db- 
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FIG. 7. The radial profile of the banana-normalized collisionality z/J and the rotational transform t. 



FIG. 8. The electric field dependence of Fi at yVVV'a — 0.29. Following abbreviations are used to represent the particle orbit 
used in a simulation; DKES-like for the particle orbit in the DKES-like limit, ZOW for the Zero Orbit Width limit, and ZMD 
for the Zero Magnetic Drift, limit. Results obtained by using DKES, GSRAKE, and FORTEC-3D codes are also plotted, where 
FORTEC-3D is denoted as Global (F3D). The ambipolar E r estimated by GSRAKE code is also shown by a vertical line. 


From Fig. 8 - 10, the particle flux of DKES-like and ZMD limits reproduce almost the same E r dependence 
as GSRAKE at every magnetic surface. Only a slight difference from original DKES code is also observed. The 
momentum conservation does not affect the local neoclassical transport due to the low collisionality of the plasma 
considered here. Also, the results of the local code in DKES-like and ZMD limits well agree with each other except for 
a very small difference in small E r of E r — 1.0 kV/m, indicating that the mono-energetic particles and incompressible 
E x B drift does not affect so much on the resultant neoclassical transport. The insignificance of these two assumptions 
are accounted for as follows. As Landreman et al. pointed out, the mono-energetic assumption varies the fraction 
of trapped- and untrapped-boundary in the velocity space, leading to the underestimation of the trapped particle 
with DKES-like orbit, 13 and the neoclassical transport in the DKES-like limit begins to give a different prediction 
from that in the ZMD limit when the poloidal mach number M p exceeds approximately 0.3. 15 Since E r considered 
here is small enough to satisfy the drift ordering of VE/vth ~ 0(<5), the difference in the trapped-particle fraction in 
DKES-like orbit (mono-energetic particles) and ZMD orbit (energy-distributed particles) does not appear so much. 

The results of ZOW limit and global code (FORTEC-3D) clearly show a different dependence of Tj on E r from that 
of other local limits and codes. The particle flux of other local simulations show a peak at E r = 0 at every surface. 
The large T; there is caused by the helically-trapped particles which have a large step size in the radial direction. 
This can be understood from the discussion around eq.(9) in Park et al. 3i Consider a case of zero tangential magnetic 
drift in the radially local system (i.e. uib = 0). All the helically-trapped particles cannot move along the surface and 
remain trapped when E T = 0 (u>e = 0), giving rise to the large radial transport. This is the reason why the radial 
flux in conventional local codes is enhanced and shows a strong peak at E T = 0. On the other hand, the particle flux 
with the ZOW orbit has no clear peak at I — 0.29 in Fig. 8, and shows small peaks at the small negative E r 

at vA'VV’a — 0.49 and 0.74 in Figs. 9 and 10, respectively. The similar tendency is also seen in the results of the 
global code. The existence of vb in ZOW limit (and the global code) makes helically-trapped particles move along the 
surface even without the E x B drift. Hence, the so-called poloidal resonance occurs at finite E r since the resonance 
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E r [kV/m] 


FIG. 9. The electric field dependence of Fi at — 0.49. An Enlarged view for —1 < E r < 1 is also shown in the upper 

right of the figure. Legends in the figure are the same as Fig. 8. The ambipolar E T estimated by GSRAKE code is also shown 
by a vertical line. 



E r [kV/m] 


FIG. 10. The electric field dependence of Ti at ^ip/if ’a. — 0.74. An Enlarged view for —1.5 < E r < 1.5 is also shown in the 
upper right of the figure. Legends in the figure are the same as Fig. 8. The ambipolar E r estimated by GSRAKE code is also 
shown by a vertical line. 


condition cob + w e = 0 is satisfied with the finite E r . Again from eq. (9) of Park et al., 3A the magnetic precession 
frequency uib becomes positive for ion species, indicating that the resonance occurs when the E x B precession is 
negative and the radial transport is enhanced at the negative E r . Moreover, the fraction of trapped particles which 
satisfies the resonance condition lob +wb = 0 is reduced due to the ^-dependence of ub, compared to the the zero 
Db case. This makes the peaked radial flux broader and smaller in the ZOW limit. 

The decrease of Tj at the resonance also arises due to the tangential drift. The complicated orbit by the combination 
of Db and E x B drifts causes a transition from trapped to untrapped particles along the surface. This is called the 
collisionless detrapping of the particle, which leads to almost no clear peak near the axis (Fig. 8), or small and broad 
peak around mid-radii (Fig. 9 and 10). Since the collisionality is low towards the plasma outer region, the effect of the 
orbit arises more significantly at ~ 0.74 than at other two surfaces. Finally, the results of FORTEC-3D show 

somewhat smaller Tj at every surface compared to those in ZOW limit, while the peak positions of these two cases are 
almost the same (E r ~ —1.6 kV/m at — 0.74). This is explained as follows. The peak position is determined 

by the balance between Db and E x B drift, and it agrees with each other since the same tangential magnetic drift 
is used in the ZOW limit and the global code. On the other hand, the finite radial drift is the only included in the 
global code. The additional effect causes another collisionless detrapping, leading to smaller Fb Also, the extent to 
which Ti shows smaller value than that of the ZOW case becomes larger towards the plasma edge in Figs. 8 - 10. This 
is attributed to the larger variation of the magnetic field experienced by the particle along the radial drift, ijidB/dip 
towards the edge in the LFID configuration, resulting in the larger fraction of the detrapped particles and smaller Ti 
of the global code than those in the ZOW limit. 
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VI. SUMMARY 

In this paper, we provide an alternative way for numerical evaluation of the local neoclassical transport with several 
kinds of the particle orbit. It aims to develop a new local neoclassical transport code which includes the finite magnetic 
drift tangential to a flux surface. Such particle orbit is called the zero orbit width (ZOW) limit in this paper since the 
radial drift is only ignored in the drift kinetic equation. The drift kinetic equation and its variations with various local 
assumptions are systematically derived from the global version with the finite orbit width (FOW) effect to the ZOW, 
zero magnetic drift (ZMD) and DKES-like limits. The systematic derivation enables us to investigate the effect of 
the tangential magnetic drift in the local neoclassical transport by comparing to the global neoclassical transport and 
other local transport models. The most significant change in the ZOW limit is that the finite tangential magnetic drift 
gives rise to the compressibility of the phase-space volume in the radially local (four-dimensional) phase space. The 
conservative property of the phase-space volume, which varies depending on the orbit and phase space considered, is 
discussed in detail. Based on the discussion of Hu and Krommes, such non-Hamiltonian (non-conservative) system 
can be appropriately treated by regarding the compressibility as a source term to the system. With this formulation, 
the two-weight 6f Monte Carlo method is presented. 

It is worth describing the difference between our source term and that proposed in Landreman et al. 15 As discussed 
in the reference, the full and partial trajectories give rise to a singular perturbation problem in surface averaged 
conservation laws of the particle number and energy when the £j- approaches to 0, where the full trajectory corresponds 
to the ZMD orbit in this paper. This is due to the fact that only the E r term survives in the conservation equations for 
the full and partial trajectory models, leading to an unphysical behavior of the distribution function in the E r = 0 limit, 
see eqs. (23) and (26) in the reference. The singular perturbation problem is successfully eliminated by introducing 
particle and/or heat source terms in the drift kinetic equation. On the other hand, in this paper, the compressibility 
of the phase-space volume arising from the finite tangential magnetic drift acts as a source term, which is of a higher 
order in the drift kinetic equation in the ZOW limit. Since the compressibility changes the conservation equations, the 
singular perturbation problem is avoided. It should be emphasized that our source term is introduced to practically 
evaluate the neoclassical transport in the ZOW model by solving the non-Hamiltonian drift kinetic equation as an 
initial value problem. Although the source term actually violates the particle number and/or energy conservations, 
in most cases presented here except for the near-axis region, it does not cause any matter in evaluating neoclassical 
particle and energy fluxes. The impact of the source term on the neoclassical transport will be discussed more in 
detail in future works. 

The code verification and its validity are checked by theoretical and numerical benchmark calculations for ax- 
isymmetric and non-axisymmetric magnetic field configurations. For a tokamak case, the neoclassical ion thermal 
diffusivity in the axisymmetric plasma well reproduces the Chang-Hinton formula in a wide range of the collisionality. 
Also, the parallel flow coefficient of the local code with DKES-like orbit and the ZOW orbit are shown to well agree 
with the theoretical estimations of Hirshman-Sigmar. This indicates that the finite magnetic drift does not change 
the conventional local neoclassical transport so much in an axisymmetric configuration. 

In a non-axisymmetric device, the finite tangential magnetic drift significantly changes the local neoclassical trans¬ 
port. In conventional local neoclassical transport calculations (ZMD and DKES-like limits), the poloidal resonance 
condition, vb — ve, where the helically-trapped particle causes a large radial transport, is satisfied with E r = 0 due to 
the absence of vb- When he exists (in the ZOW limit), however, the poloidal resonance is shifted to a small negative 
E r . It is demonstrated for the first time that the finite tangential magnetic drift gives rise to a qualitative change in 
the dependence of the neoclassical transport on the radial electric field even in the local neoclassical transport model. 
Also, similarly to the global neoclassical transport, the large radial transport at the resonance seen in the ZMD and 
DKES-like limits can be avoided in the ZOW limit due to the collision detrapping along the local orbit. Hence, two 
important physics included in the global code can be captured by the local code in the ZOW limit. 

A key role of the neoclassical transport in a non-axisymmetric plasma is to predict the ambipolar K r according 
to the ambipolar condition of the neoclassical particle flux. As demonstrated in the paper, the finite magnetic drift 
changes the dependence of the ion particle flux on E r . The ambipolar E t predicted can vary depending on whether 
Vb is included in evaluating the neoclassical transport. The main cause of the difference comes from the shift of the 
poloidal resonance condition, and it is included in our local code in the ZOW limit. Since the local code is less time- 
consuming and requires less computational cost than the global code, the local code will be a preferable alternative 
to predict the ambipolar E z in experimental analyses. 
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